Analysis of simulated vs observed runoff
#Graph observed vs simulated runoff for Max Multiplicative Change Factor (MCF)
observed<- read.csv(here("climate_change_data","SWMM_results", "Mar14_observed.csv"))
simulated_max_march14 <- read.csv(here("climate_change_data","SWMM_results", "MCF_max_Mar14_simulated.csv"))
simulated_min_march14 <- read.csv(here("climate_change_data","SWMM_results", "MCF_min_Mar14_simulated.csv"))
# merge MCF max with observed
discharge_MCF_max<- merge(observed, simulated_max_march14, by = "time_step")
#merge MCF min with observed
discharge_MCF_min<- merge(observed, simulated_min_march14, by = "time_step")
summary(discharge_MCF_max)
# statical analyses
# dischargettest<- t.test(discharge_MCF_max$discharge_obs_cfs, discharge_MCF_max$simulated_flow_cfs)
#
# regress<- lm(discharge_obs_cfs ~ simulated_flow_cfs, data = discharge_MCF_max)
#
# regress
#
# sim<- discharge_MCF_max$simulated_flow_cfs
# obs<- discharge_MCF_max$discharge_obs_cfs
#
#
# sutcliffe<- NSE(sim, obs, na.rm=TRUE, FUN=NULL, epsilon=c("Pushpalatha2012"))
#
# sutcliffe
# write.csv(discharge,file = "discharge_obv_sim1.csv", row.names = FALSE)
graph_max<- ggplot(discharge_MCF_max, aes(x=discharge_obs_cfs, simulated_flow_cfs))+
geom_point()
graph_max
graph_min<- ggplot(discharge_MCF_min, aes(x=discharge_obs_cfs, simulated_flow_cfs))+
geom_point()
graph_min
calibrate_graph_max<- discharge_MCF_max %>%
ggplot()+
geom_line(aes(x=time_step, y=discharge_obs_cfs), color="#000000", size = 2.5)+
geom_line(aes(x=time_step, y=simulated_flow_cfs), color= "#009E73", size = 2.5)+
theme_classic()+
labs(x="Time (hours)", y="Discharge (cfs)")+
scale_y_continuous(limits= c(0,850), breaks= seq(0,850, by= 100),expand= c(0,0))+
scale_x_continuous(limits= c(0,12), breaks= seq(0,12, by= 2),expand= c(0,0))+
annotate("text", label= "Observed", x=10.5, y=85, size=8)+
annotate("text", label= "Simulated", x=10.5, y=300, size=8) +
theme(text = element_text(size=22))
calibrate_graph_max
calibrate_graph_min<- discharge_MCF_min %>%
ggplot()+
geom_line(aes(x=time_step, y=discharge_obs_cfs), color="#000000", size =2.5) +
geom_line(aes(x=time_step, y=simulated_flow_cfs), color= "#009E73", size =2.5)+
theme_classic()+
labs(x="Time (hours)", y="Discharge (cfs)")+
scale_y_continuous(limits= c(0,400), breaks= seq(0,350, by= 50),expand= c(0,0))+
scale_x_continuous(limits= c(0,12), breaks= seq(0,12, by= 2),expand= c(0,0))+
annotate("text", label= "Simulated", x=10.5, y=45, size=8)+
annotate("text", label= "Observed", x=10.5, y=185, size=8) +
theme(text = element_text(size=22))
calibrate_graph_min
# save graphs
# ggsave('discharge_mar14_max_MCF.png', calibrate_graph_max, width = 16, height = 9, units = "in")
# ggsave('discharge_mar14_min_MCF.png', calibrate_graph_min, width = 16, height = 9, units = "in")
## combine all three lines onto one graph for easabilty
calibrate_graph <- ggplot()+
geom_line(data = discharge_MCF_min, aes(x=time_step, y=discharge_obs_cfs), color="#000000", size =2.5) +
geom_line(data= discharge_MCF_min, aes(x=time_step, y=simulated_flow_cfs), color= "#009E73", size =2.5)+
geom_line(data= discharge_MCF_max, aes(x=time_step, y=simulated_flow_cfs), color= "#D55E00", size = 2.5) +
theme_classic()+
labs(x="Time (hours)", y="Discharge (cfs)")+
scale_y_continuous(limits= c(0,850), breaks= seq(0,850, by= 100),expand= c(0,0))+
scale_x_continuous(limits= c(0,12), breaks= seq(0,12, by= 2),expand= c(0,0))+
annotate("text", label= "Observed", x=10.5, y=200, size=8)+
annotate("text", label= "Max MCF Simulated", x=10.5, y=655, size=8) +
annotate("text", label= "Min MCF Simulated", x=10.5, y=55, size=8) +
theme(text = element_text(size=22))
calibrate_graph
# ggsave('discharge_mar14_MCF.png', calibrate_graph, width = 16, height = 9, units = "in")
Analysis of MCF max storm results
subcatch_max_mcf<- read_csv("subcatchments_all.csv") %>%
mutate(subcatchment= OBJECTID_1) %>%
select(subcatchment, Curve_Number, Slope, percent_imp, Area_sqft)
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## OBJECTID_1 = col_double(),
## Area_sqft = col_double(),
## Curve_Number = col_double(),
## Slope = col_double(),
## percent_imp = col_double()
## )
swm_results_max_mcf<-read_csv("Wailupe_MCF_Max_Mar14.csv")
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## subcatchment = col_double(),
## total_precip_in = col_double(),
## total_runon_in = col_double(),
## total_evap_in = col_double(),
## total_infil_in = col_double(),
## imperv_runoff_in = col_double(),
## perv_runoff_in = col_double(),
## total_runoff_in = col_double(),
## total_runoff_10_6_gal = col_double(),
## peak_runoff_cfs = col_double(),
## runoff_coeff = col_double()
## )
#Characterize results by urbanization level
results_max_mcf<- merge(subcatch_max_mcf, swm_results_max_mcf, by = "subcatchment") %>%
mutate(runoff_normalized=
total_runoff_in/Area_sqft) %>%
mutate(Urbanization_level=
case_when(
percent_imp <15 |percent_imp == 15 ~ "Natural (less than 15 % Impervious)",
percent_imp >15 & percent_imp <45 ~ "Urbanized (Between 15 and 45 % Impervious)",
percent_imp >44.9999 ~ "Very urbanized (More than 45 % Impervious)"
)
)
write.csv(results_max_mcf, file = "results_max_mcf.csv") ##Export as .csv for use with graph maps
## normalize peak discharge (cfs) for max MCF by dividing by subcatchment area
results_max_mcf_normalized <- results_max_mcf %>%
mutate(peak_runoff_cfs_norm = peak_runoff_cfs/Area_sqft)
results_maxmcf_peakflow_table <- results_max_mcf %>%
select(subcatchment, peak_runoff_cfs, Curve_Number, Slope, percent_imp, Area_sqft) %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
arrange(desc(peak_runoff_cfs)) %>%
top_n(20,peak_runoff_cfs) %>%
kable(col.names=c("Subcatchment","Peak Runoff (cfs), Max MCF",
"Curve Number", "Slope",
"Percent Impervious", "Area (sqft)")) %>%
kable_styling(bootstrap_options = "striped")
results_maxmcf_peakflow_table
|
Subcatchment
|
Peak Runoff (cfs), Max MCF
|
Curve Number
|
Slope
|
Percent Impervious
|
Area (sqft)
|
|
63
|
106.39
|
55.59826
|
11.3508475
|
59.595082
|
1260906.0
|
|
11
|
105.54
|
55.59506
|
0.9475857
|
60.176841
|
1461121.3
|
|
28
|
99.65
|
58.22901
|
9.6951530
|
45.648739
|
1400355.0
|
|
74
|
90.09
|
45.63212
|
0.5130342
|
49.951879
|
1593473.1
|
|
42
|
80.56
|
49.94638
|
12.7741916
|
52.075244
|
1097149.0
|
|
89
|
79.13
|
58.82177
|
12.2672167
|
54.379880
|
929708.8
|
|
78
|
69.55
|
42.51236
|
0.6625196
|
45.869531
|
1272102.3
|
|
20
|
64.99
|
68.12691
|
28.7382083
|
5.926749
|
3260351.4
|
|
62
|
62.88
|
47.11242
|
4.9693069
|
57.621143
|
816823.7
|
|
23
|
59.82
|
58.75298
|
2.0917640
|
60.158489
|
712111.1
|
|
58
|
57.25
|
46.82201
|
13.4737125
|
50.489651
|
794612.3
|
|
7
|
56.98
|
62.46735
|
24.3708501
|
14.101788
|
2174275.2
|
|
76
|
52.96
|
61.94713
|
12.6689327
|
21.946154
|
1533842.8
|
|
79
|
48.58
|
42.64164
|
0.7809627
|
55.674414
|
703208.6
|
|
75
|
47.57
|
57.16000
|
3.1037990
|
49.969260
|
646779.6
|
|
91
|
47.33
|
61.92057
|
7.2257735
|
49.197766
|
543493.1
|
|
10
|
38.42
|
59.50282
|
13.0349991
|
46.827316
|
527259.5
|
|
52
|
38.27
|
46.04250
|
14.2092259
|
49.068201
|
521285.7
|
|
67
|
37.74
|
59.21465
|
13.0945141
|
65.315445
|
376961.0
|
|
34
|
35.34
|
50.45073
|
6.5716791
|
55.845327
|
439337.8
|
results_maxmcf_peakflow_table_normalized <- results_max_mcf_normalized %>%
select(subcatchment, peak_runoff_cfs, Curve_Number, Slope, percent_imp, Area_sqft,peak_runoff_cfs_norm ) %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
arrange(desc(peak_runoff_cfs_norm)) %>%
top_n(20,peak_runoff_cfs_norm) %>%
kable(col.names=c("Subcatchment","Peak Runoff (cfs), Max MCF",
"Curve Number", "Slope",
"Percent Impervious", "Area (sqft)", "Normalized Peak Runoff (cfs/sqft)")) %>%
kable_styling(bootstrap_options = "striped")
results_maxmcf_peakflow_table_normalized
|
Subcatchment
|
Peak Runoff (cfs), Max MCF
|
Curve Number
|
Slope
|
Percent Impervious
|
Area (sqft)
|
Normalized Peak Runoff (cfs/sqft)
|
|
5
|
4.51
|
59.94459
|
4.973908
|
59.22529
|
42008.07
|
0.0001074
|
|
60
|
13.90
|
58.26651
|
10.200056
|
66.27512
|
130703.86
|
0.0001063
|
|
68
|
31.68
|
64.12348
|
16.778001
|
56.40240
|
306366.86
|
0.0001034
|
|
45
|
21.10
|
46.68192
|
10.837372
|
75.81188
|
205273.48
|
0.0001028
|
|
46
|
21.33
|
57.85258
|
9.180995
|
59.97949
|
209977.29
|
0.0001016
|
|
47
|
16.51
|
44.12992
|
9.017865
|
73.67874
|
163180.85
|
0.0001012
|
|
67
|
37.74
|
59.21465
|
13.094514
|
65.31544
|
376961.02
|
0.0001001
|
|
40
|
7.68
|
46.96258
|
7.774446
|
65.26444
|
78272.63
|
0.0000981
|
|
71
|
15.57
|
58.22480
|
1.229756
|
63.33945
|
161251.17
|
0.0000966
|
|
14
|
23.52
|
59.93120
|
10.062202
|
51.97367
|
244581.39
|
0.0000962
|
|
22
|
10.23
|
52.60335
|
5.224538
|
57.36459
|
106889.43
|
0.0000957
|
|
69
|
16.14
|
61.92279
|
23.710094
|
39.74707
|
171176.11
|
0.0000943
|
|
65
|
29.11
|
59.79150
|
12.285943
|
54.98557
|
309782.87
|
0.0000940
|
|
29
|
5.84
|
46.60531
|
4.303893
|
63.75090
|
62214.52
|
0.0000939
|
|
31
|
21.49
|
58.70699
|
10.270774
|
47.00970
|
229731.45
|
0.0000935
|
|
59
|
22.41
|
53.74048
|
9.355130
|
58.61615
|
241399.55
|
0.0000928
|
|
82
|
5.99
|
52.36590
|
12.984541
|
50.03241
|
65510.58
|
0.0000914
|
|
51
|
17.68
|
45.48179
|
3.205696
|
65.70618
|
194382.01
|
0.0000910
|
|
54
|
27.29
|
46.30916
|
11.985041
|
65.90850
|
301206.44
|
0.0000906
|
|
21
|
7.00
|
41.71682
|
2.133592
|
64.38962
|
78832.74
|
0.0000888
|
results_maxmcf_runoffcoef_table <- results_max_mcf %>%
select(subcatchment, runoff_coeff, Curve_Number, Slope, percent_imp, Area_sqft) %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
arrange(desc(runoff_coeff)) %>%
top_n(20,runoff_coeff) %>%
kable(col.names=c("Subcatchment","Runoff Coefficient, Max MCF",
"Curve Number", "Slope",
"Percent Impervious", "Area (sqft)")) %>%
kable_styling(bootstrap_options = "striped")
results_maxmcf_runoffcoef_table
|
Subcatchment
|
Runoff Coefficient, Max MCF
|
Curve Number
|
Slope
|
Percent Impervious
|
Area (sqft)
|
|
60
|
0.821
|
58.26651
|
10.200056
|
66.27512
|
130703.86
|
|
45
|
0.815
|
46.68192
|
10.837372
|
75.81188
|
205273.48
|
|
68
|
0.813
|
64.12348
|
16.778001
|
56.40240
|
306366.86
|
|
5
|
0.809
|
59.94459
|
4.973908
|
59.22529
|
42008.07
|
|
67
|
0.809
|
59.21465
|
13.094514
|
65.31544
|
376961.02
|
|
47
|
0.799
|
44.12992
|
9.017865
|
73.67874
|
163180.85
|
|
71
|
0.793
|
58.22480
|
1.229756
|
63.33945
|
161251.17
|
|
46
|
0.789
|
57.85258
|
9.180995
|
59.97949
|
209977.29
|
|
65
|
0.770
|
59.79150
|
12.285943
|
54.98557
|
309782.87
|
|
14
|
0.764
|
59.93120
|
10.062202
|
51.97367
|
244581.39
|
|
40
|
0.762
|
46.96258
|
7.774446
|
65.26444
|
78272.63
|
|
23
|
0.759
|
58.75298
|
2.091764
|
60.15849
|
712111.05
|
|
91
|
0.755
|
61.92057
|
7.225773
|
49.19777
|
543493.09
|
|
59
|
0.751
|
53.74048
|
9.355130
|
58.61615
|
241399.55
|
|
54
|
0.748
|
46.30916
|
11.985041
|
65.90850
|
301206.44
|
|
89
|
0.748
|
58.82177
|
12.267217
|
54.37988
|
929708.78
|
|
29
|
0.747
|
46.60531
|
4.303893
|
63.75090
|
62214.52
|
|
51
|
0.746
|
45.48179
|
3.205696
|
65.70618
|
194382.01
|
|
63
|
0.746
|
55.59826
|
11.350848
|
59.59508
|
1260905.98
|
|
22
|
0.745
|
52.60335
|
5.224538
|
57.36459
|
106889.43
|
#Perform a linear regression for the SWMM max MCF storm results
results_max_mcf_regression<- lm(total_runoff_in ~Curve_Number + Slope + percent_imp +
Area_sqft, data = results_max_mcf)
#Graph the relationship bw simulated runoff and impervious cover by urbanization level
runoff_imp_graph_max_mcf<- results_max_mcf %>%
ggplot(aes(x=percent_imp, y=total_runoff_in))+
geom_point(aes(color=Urbanization_level))+
labs(x= "Percent Impervious of Subcatchment", y= "Total Simulated Runoff (inches)")+
scale_y_continuous(limits= c(0,8), breaks= seq(0,8, by= 1),expand= c(0,0.08))+
scale_x_continuous(limits= c(0,80), breaks= seq(0,80, by= 10),expand= c(0,0))+
scale_color_manual(name= "Urbanization Level", values= c("darkgreen", "darkseagreen",
"darkgoldenrod1"))+
theme_classic()
runoff_imp_graph_max_mcf
#Save runoff vs. impervious graph
# ggsave("runoff_imp_max_mcf.pdf", width = 8, height =4)
# ggsave("runoff_imp_max_mcf.png", width = 8, height =4)
#Create a table for the regression results for max mcf storm
regress_table_max_mcf<- stargazer(results_max_mcf_regression, type ="html", digits= 2,
dep.var.labels = "Total Runoff (Inches)",
covariate.labels = c("Curve Number", "Slope", "Percent Impervious",
"Area (sqft)", "Y-Intercept"),
omit.stat = c("rsq"))
##
## <table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>Total Runoff (Inches)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Curve Number</td><td>0.03<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.01)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Slope</td><td>0.05<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.01)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Percent Impervious</td><td>0.05<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.01)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Area (sqft)</td><td>0.0000<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.0000)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Y-Intercept</td><td>-0.34</td></tr>
## <tr><td style="text-align:left"></td><td>(0.61)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>97</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.48</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>0.67 (df = 92)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>23.33<sup>***</sup> (df = 4; 92)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
regress_table_max_mcf
## [1] ""
## [2] "<table style=\"text-align:center\"><tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\"></td><td><em>Dependent variable:</em></td></tr>"
## [3] "<tr><td></td><td colspan=\"1\" style=\"border-bottom: 1px solid black\"></td></tr>"
## [4] "<tr><td style=\"text-align:left\"></td><td>Total Runoff (Inches)</td></tr>"
## [5] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\">Curve Number</td><td>0.03<sup>**</sup></td></tr>"
## [6] "<tr><td style=\"text-align:left\"></td><td>(0.01)</td></tr>"
## [7] "<tr><td style=\"text-align:left\"></td><td></td></tr>"
## [8] "<tr><td style=\"text-align:left\">Slope</td><td>0.05<sup>***</sup></td></tr>"
## [9] "<tr><td style=\"text-align:left\"></td><td>(0.01)</td></tr>"
## [10] "<tr><td style=\"text-align:left\"></td><td></td></tr>"
## [11] "<tr><td style=\"text-align:left\">Percent Impervious</td><td>0.05<sup>***</sup></td></tr>"
## [12] "<tr><td style=\"text-align:left\"></td><td>(0.01)</td></tr>"
## [13] "<tr><td style=\"text-align:left\"></td><td></td></tr>"
## [14] "<tr><td style=\"text-align:left\">Area (sqft)</td><td>0.0000<sup>***</sup></td></tr>"
## [15] "<tr><td style=\"text-align:left\"></td><td>(0.0000)</td></tr>"
## [16] "<tr><td style=\"text-align:left\"></td><td></td></tr>"
## [17] "<tr><td style=\"text-align:left\">Y-Intercept</td><td>-0.34</td></tr>"
## [18] "<tr><td style=\"text-align:left\"></td><td>(0.61)</td></tr>"
## [19] "<tr><td style=\"text-align:left\"></td><td></td></tr>"
## [20] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\">Observations</td><td>97</td></tr>"
## [21] "<tr><td style=\"text-align:left\">Adjusted R<sup>2</sup></td><td>0.48</td></tr>"
## [22] "<tr><td style=\"text-align:left\">Residual Std. Error</td><td>0.67 (df = 92)</td></tr>"
## [23] "<tr><td style=\"text-align:left\">F Statistic</td><td>23.33<sup>***</sup> (df = 4; 92)</td></tr>"
## [24] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\"><em>Note:</em></td><td style=\"text-align:right\"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>"
## [25] "</table>"
Analysis of MCF min storm results
subcatch_min_mcf<- read_csv("subcatchments_all.csv") %>%
mutate(subcatchment= OBJECTID_1) %>%
select(subcatchment, Curve_Number, Slope, percent_imp, Area_sqft)
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## OBJECTID_1 = col_double(),
## Area_sqft = col_double(),
## Curve_Number = col_double(),
## Slope = col_double(),
## percent_imp = col_double()
## )
swm_results_min_mcf<-read_csv("Wailupe_MCF_Min_Mar14.csv")
##
## ── Column specification ────────────────────────────────────────────────────────
## cols(
## subcatchment = col_double(),
## total_precip_in = col_double(),
## total_runon_in = col_double(),
## total_evap_in = col_double(),
## total_infil_in = col_double(),
## imperv_runoff_in = col_double(),
## perv_runoff_in = col_double(),
## total_runoff_in = col_double(),
## total_runoff_10_6_gal = col_double(),
## peak_runoff_cfs = col_double(),
## runoff_coeff = col_double()
## )
#Characterize results by urbanization level
results_min_mcf<- merge(subcatch_min_mcf, swm_results_min_mcf, by = "subcatchment") %>%
mutate(runoff_normalized=
total_runoff_in/Area_sqft) %>%
mutate(Urbanization_level=
case_when(
percent_imp <15 |percent_imp == 15 ~ "Natural (less than 15 % Impervious)",
percent_imp >15 & percent_imp <45 ~ "Urbanized (Between 15 and 45 % Impervious)",
percent_imp >44.9999 ~ "Very urbanized (More than 45 % Impervious)"
)
)
write.csv(results_min_mcf, file = "results_min_mcf.csv") ##Export as .csv for use with graph maps
## normalize peak discharge (cfs) for min MCF by dividing by subcatchment area
results_min_mcf_normalized <- results_min_mcf %>%
mutate(peak_runoff_cfs_norm = peak_runoff_cfs/Area_sqft)
results_minmcf_peakflow_table <- results_min_mcf %>%
select(subcatchment, peak_runoff_cfs, Curve_Number, Slope, percent_imp, Area_sqft) %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
arrange(desc(peak_runoff_cfs)) %>%
top_n(20,peak_runoff_cfs) %>%
kable(col.names=c("Subcatchment","Peak Runoff (cfs), Min MCF",
"Curve Number", "Slope",
"Percent Impervious", "Area (sqft)")) %>%
kable_styling(bootstrap_options = "striped")
results_minmcf_peakflow_table
|
Subcatchment
|
Peak Runoff (cfs), Min MCF
|
Curve Number
|
Slope
|
Percent Impervious
|
Area (sqft)
|
|
63
|
15.29
|
55.59826
|
11.3508475
|
59.59508
|
1260906.0
|
|
28
|
13.56
|
58.22901
|
9.6951530
|
45.64874
|
1400355.0
|
|
42
|
12.15
|
49.94638
|
12.7741916
|
52.07524
|
1097149.0
|
|
11
|
11.91
|
55.59506
|
0.9475857
|
60.17684
|
1461121.3
|
|
89
|
10.85
|
58.82177
|
12.2672167
|
54.37988
|
929708.8
|
|
74
|
10.07
|
45.63212
|
0.5130342
|
49.95188
|
1593473.1
|
|
62
|
9.66
|
47.11242
|
4.9693069
|
57.62114
|
816823.7
|
|
58
|
8.76
|
46.82201
|
13.4737125
|
50.48965
|
794612.3
|
|
78
|
8.38
|
42.51236
|
0.6625196
|
45.86953
|
1272102.3
|
|
23
|
8.09
|
58.75298
|
2.0917640
|
60.15849
|
712111.1
|
|
76
|
7.40
|
61.94713
|
12.6689327
|
21.94615
|
1533842.8
|
|
7
|
6.79
|
62.46735
|
24.3708501
|
14.10179
|
2174275.2
|
|
75
|
6.65
|
57.16000
|
3.1037990
|
49.96926
|
646779.6
|
|
79
|
6.63
|
42.64164
|
0.7809627
|
55.67441
|
703208.6
|
|
91
|
5.86
|
61.92057
|
7.2257735
|
49.19777
|
543493.1
|
|
52
|
5.66
|
46.04250
|
14.2092259
|
49.06820
|
521285.7
|
|
10
|
5.46
|
59.50282
|
13.0349991
|
46.82732
|
527259.5
|
|
67
|
5.35
|
59.21465
|
13.0945141
|
65.31544
|
376961.0
|
|
34
|
5.31
|
50.45073
|
6.5716791
|
55.84533
|
439337.8
|
|
55
|
4.37
|
52.00222
|
14.0124563
|
43.23872
|
456026.8
|
results_minmcf_peakflow_normalized_table <- results_min_mcf_normalized %>%
select(subcatchment, peak_runoff_cfs, Curve_Number, Slope, percent_imp, Area_sqft,peak_runoff_cfs_norm ) %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
arrange(desc(peak_runoff_cfs_norm)) %>%
top_n(20,peak_runoff_cfs_norm) %>%
kable(col.names=c("Subcatchment","Peak Runoff (cfs), Min MCF",
"Curve Number", "Slope",
"Percent Impervious", "Area (sqft)", "Normalized Peak Runoff (cfs/sqft)")) %>%
kable_styling(bootstrap_options = "striped")
results_minmcf_peakflow_normalized_table
|
Subcatchment
|
Peak Runoff (cfs), Min MCF
|
Curve Number
|
Slope
|
Percent Impervious
|
Area (sqft)
|
Normalized Peak Runoff (cfs/sqft)
|
|
45
|
3.41
|
46.68192
|
10.837372
|
75.81188
|
205273.48
|
0.0000166
|
|
47
|
2.66
|
44.12992
|
9.017865
|
73.67874
|
163180.85
|
0.0000163
|
|
60
|
1.92
|
58.26651
|
10.200056
|
66.27512
|
130703.86
|
0.0000147
|
|
40
|
1.13
|
46.96258
|
7.774446
|
65.26444
|
78272.63
|
0.0000144
|
|
51
|
2.80
|
45.48179
|
3.205696
|
65.70618
|
194382.01
|
0.0000144
|
|
54
|
4.32
|
46.30916
|
11.985041
|
65.90850
|
301206.44
|
0.0000143
|
|
21
|
1.12
|
41.71682
|
2.133592
|
64.38962
|
78832.74
|
0.0000142
|
|
67
|
5.35
|
59.21465
|
13.094514
|
65.31544
|
376961.02
|
0.0000142
|
|
29
|
0.88
|
46.60531
|
4.303893
|
63.75090
|
62214.52
|
0.0000141
|
|
71
|
2.21
|
58.22480
|
1.229756
|
63.33945
|
161251.17
|
0.0000137
|
|
49
|
4.29
|
44.95489
|
11.626618
|
61.33458
|
317060.05
|
0.0000135
|
|
46
|
2.80
|
57.85258
|
9.180995
|
59.97949
|
209977.29
|
0.0000133
|
|
13
|
2.77
|
37.24907
|
2.801558
|
60.72206
|
210147.80
|
0.0000132
|
|
5
|
0.55
|
59.94459
|
4.973908
|
59.22529
|
42008.07
|
0.0000131
|
|
59
|
3.13
|
53.74048
|
9.355130
|
58.61615
|
241399.55
|
0.0000130
|
|
35
|
2.11
|
36.38429
|
2.502320
|
58.10591
|
163977.29
|
0.0000129
|
|
22
|
1.36
|
52.60335
|
5.224538
|
57.36459
|
106889.43
|
0.0000127
|
|
36
|
1.31
|
36.48504
|
3.413503
|
57.37835
|
103198.24
|
0.0000127
|
|
44
|
3.09
|
46.81877
|
9.792575
|
56.58471
|
246676.73
|
0.0000125
|
|
68
|
3.83
|
64.12348
|
16.778001
|
56.40240
|
306366.86
|
0.0000125
|
results_minmcf_runoffcoef_table <- results_min_mcf %>%
select(subcatchment, runoff_coeff, Curve_Number, Slope, percent_imp, Area_sqft) %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
arrange(desc(runoff_coeff)) %>%
top_n(20,runoff_coeff) %>%
kable(col.names=c("Subcatchment","Runoff Coefficient, Min MCF",
"Curve Number", "Slope",
"Percent Impervious", "Area (sqft)")) %>%
kable_styling(bootstrap_options = "striped")
results_minmcf_runoffcoef_table
|
Subcatchment
|
Runoff Coefficient, Min MCF
|
Curve Number
|
Slope
|
Percent Impervious
|
Area (sqft)
|
|
45
|
0.594
|
46.68192
|
10.837372
|
75.81188
|
205273.48
|
|
47
|
0.586
|
44.12992
|
9.017865
|
73.67874
|
163180.85
|
|
60
|
0.531
|
58.26651
|
10.200056
|
66.27512
|
130703.86
|
|
40
|
0.528
|
46.96258
|
7.774446
|
65.26444
|
78272.63
|
|
51
|
0.516
|
45.48179
|
3.205696
|
65.70618
|
194382.01
|
|
54
|
0.514
|
46.30916
|
11.985041
|
65.90850
|
301206.44
|
|
21
|
0.511
|
41.71682
|
2.133592
|
64.38962
|
78832.74
|
|
29
|
0.510
|
46.60531
|
4.303893
|
63.75090
|
62214.52
|
|
67
|
0.509
|
59.21465
|
13.094514
|
65.31544
|
376961.02
|
|
71
|
0.493
|
58.22480
|
1.229756
|
63.33945
|
161251.17
|
|
5
|
0.489
|
59.94459
|
4.973908
|
59.22529
|
42008.07
|
|
49
|
0.484
|
44.95489
|
11.626618
|
61.33458
|
317060.05
|
|
46
|
0.481
|
57.85258
|
9.180995
|
59.97949
|
209977.29
|
|
13
|
0.473
|
37.24907
|
2.801558
|
60.72206
|
210147.80
|
|
59
|
0.465
|
53.74048
|
9.355130
|
58.61615
|
241399.55
|
|
22
|
0.463
|
52.60335
|
5.224538
|
57.36459
|
106889.43
|
|
35
|
0.461
|
36.38429
|
2.502320
|
58.10591
|
163977.29
|
|
36
|
0.458
|
36.48504
|
3.413503
|
57.37835
|
103198.24
|
|
63
|
0.454
|
55.59826
|
11.350848
|
59.59508
|
1260905.98
|
|
23
|
0.452
|
58.75298
|
2.091764
|
60.15849
|
712111.05
|
#Perform a linear regression for the SWMM max MCF storm results
results_min_mcf_regression<- lm(total_runoff_in ~Curve_Number + Slope + percent_imp +
Area_sqft, data = results_min_mcf)
#Graph the relationship bw simulated runoff and impervious cover by urbanization level
runoff_imp_graph_min_mcf<- results_min_mcf %>%
ggplot(aes(x=percent_imp, y=total_runoff_in))+
geom_point(aes(color=Urbanization_level))+
labs(x= "Percent Impervious of Subcatchment", y= "Total Simulated Runoff (inches)")+
scale_y_continuous(limits= c(0,7), breaks= seq(0,7, by= 1),expand= c(0,0.08))+
scale_x_continuous(limits= c(0,80), breaks= seq(0,80, by= 10),expand= c(0,0))+
scale_color_manual(name= "Urbanization Level", values= c("darkgreen", "darkseagreen",
"darkgoldenrod1"))+
theme_classic()
runoff_imp_graph_min_mcf
#Save runoff vs. impervious graph
# ggsave("runoff_imp_min_mcf.pdf", width = 8, height =4)
# ggsave("runoff_imp_min_mcf.png", width = 8, height =4)
#Create a table for the regression results
regress_table_min_mcf<- stargazer(results_min_mcf_regression, type ="html", digits= 2,
dep.var.labels = "Total Runoff (Inches)",
covariate.labels = c("Curve Number", "Slope", "Percent Impervious",
"Area (sqft)", "Y-Intercept"),
omit.stat = c("rsq"))
##
## <table style="text-align:center"><tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"></td><td><em>Dependent variable:</em></td></tr>
## <tr><td></td><td colspan="1" style="border-bottom: 1px solid black"></td></tr>
## <tr><td style="text-align:left"></td><td>Total Runoff (Inches)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Curve Number</td><td>0.001</td></tr>
## <tr><td style="text-align:left"></td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Slope</td><td>0.002<sup>*</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.001)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Percent Impervious</td><td>0.01<sup>***</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.0005)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Area (sqft)</td><td>0.00<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.00)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td style="text-align:left">Y-Intercept</td><td>-0.13<sup>**</sup></td></tr>
## <tr><td style="text-align:left"></td><td>(0.05)</td></tr>
## <tr><td style="text-align:left"></td><td></td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left">Observations</td><td>97</td></tr>
## <tr><td style="text-align:left">Adjusted R<sup>2</sup></td><td>0.90</td></tr>
## <tr><td style="text-align:left">Residual Std. Error</td><td>0.06 (df = 92)</td></tr>
## <tr><td style="text-align:left">F Statistic</td><td>220.89<sup>***</sup> (df = 4; 92)</td></tr>
## <tr><td colspan="2" style="border-bottom: 1px solid black"></td></tr><tr><td style="text-align:left"><em>Note:</em></td><td style="text-align:right"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>
## </table>
regress_table_min_mcf
## [1] ""
## [2] "<table style=\"text-align:center\"><tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\"></td><td><em>Dependent variable:</em></td></tr>"
## [3] "<tr><td></td><td colspan=\"1\" style=\"border-bottom: 1px solid black\"></td></tr>"
## [4] "<tr><td style=\"text-align:left\"></td><td>Total Runoff (Inches)</td></tr>"
## [5] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\">Curve Number</td><td>0.001</td></tr>"
## [6] "<tr><td style=\"text-align:left\"></td><td>(0.001)</td></tr>"
## [7] "<tr><td style=\"text-align:left\"></td><td></td></tr>"
## [8] "<tr><td style=\"text-align:left\">Slope</td><td>0.002<sup>*</sup></td></tr>"
## [9] "<tr><td style=\"text-align:left\"></td><td>(0.001)</td></tr>"
## [10] "<tr><td style=\"text-align:left\"></td><td></td></tr>"
## [11] "<tr><td style=\"text-align:left\">Percent Impervious</td><td>0.01<sup>***</sup></td></tr>"
## [12] "<tr><td style=\"text-align:left\"></td><td>(0.0005)</td></tr>"
## [13] "<tr><td style=\"text-align:left\"></td><td></td></tr>"
## [14] "<tr><td style=\"text-align:left\">Area (sqft)</td><td>0.00<sup>**</sup></td></tr>"
## [15] "<tr><td style=\"text-align:left\"></td><td>(0.00)</td></tr>"
## [16] "<tr><td style=\"text-align:left\"></td><td></td></tr>"
## [17] "<tr><td style=\"text-align:left\">Y-Intercept</td><td>-0.13<sup>**</sup></td></tr>"
## [18] "<tr><td style=\"text-align:left\"></td><td>(0.05)</td></tr>"
## [19] "<tr><td style=\"text-align:left\"></td><td></td></tr>"
## [20] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\">Observations</td><td>97</td></tr>"
## [21] "<tr><td style=\"text-align:left\">Adjusted R<sup>2</sup></td><td>0.90</td></tr>"
## [22] "<tr><td style=\"text-align:left\">Residual Std. Error</td><td>0.06 (df = 92)</td></tr>"
## [23] "<tr><td style=\"text-align:left\">F Statistic</td><td>220.89<sup>***</sup> (df = 4; 92)</td></tr>"
## [24] "<tr><td colspan=\"2\" style=\"border-bottom: 1px solid black\"></td></tr><tr><td style=\"text-align:left\"><em>Note:</em></td><td style=\"text-align:right\"><sup>*</sup>p<0.1; <sup>**</sup>p<0.05; <sup>***</sup>p<0.01</td></tr>"
## [25] "</table>"
Maps of SWMM Results
Code setup - Load packages
library(sf) #For shapefiles
library(tmap) #For mapmaking
library(tmaptools) #For mapmaking
library(janitor) #For cleaning names
Maps for MCF max storm hotspots
results_max_mcf<- read_csv("results_max_mcf.csv") ##read in file from code above
#Combine subcatchments outline with wet storm results
subcatch_max <- read_sf(dsn = here("climate_change_data", "SWMM_results","shapefiles"), layer = "subcatch_outline") %>%
st_transform(st_crs(4326)) %>%
clean_names() %>%
select(subcatchment = objectid_1) %>%
merge(results_max_mcf) %>%
filter(subcatchment != "5")
# Total volume hotspots
hotspots_max_mcf_total <- tm_basemap("OpenStreetMap.Mapnik") +
tm_shape(subcatch_max, unit = "Miles") +
tm_polygons("runoff_coeff", alpha = 0.8, palette = "Blues", style = "cont", n=8,
legend.hist = TRUE, title = "Runoff Coefficient") +
tm_layout(title = "March 2009 Max MCF storm", inner.margins=c(.05, .05, 0.1, .53),
legend.position = c(.6,.32), legend.title.size = 1.4, legend.text.size = 1) +
tm_text("subcatchment", size = 0.3) +
tm_scale_bar(position = c(.6,.59), breaks = c(0, 0.2, 0.4, 0.6, 0.8,1)) +
tm_compass(position = c(.58,.52))
# tmap_save(hotspots_max_mcf_total, here("climate_change_data", "output_maps","hotspots_mcf_max.png")) ## not exporting
# Peak flow hotspots
hotspots_max_mcf_peak <- tm_basemap("OpenStreetMap.Mapnik") +
tm_shape(subcatch_max, unit = "Miles") +
tm_polygons("peak_runoff_cfs", alpha = 0.75, palette = "Greens", style = "cont", n=8,
legend.hist = TRUE, title = "Peak Discharge (cfs)") +
tm_layout(title = "March 2009 Max MCF storm", inner.margins=c(.05, .05, 0.1, .53),
legend.position = c(.6,.27), legend.title.size = 1.4, legend.text.size = 1) +
tm_text("subcatchment", size = 0.3) +
tm_scale_bar(position = c(.6,.54), breaks = c(0, 0.2, 0.4, 0.6, 0.8,1)) +
tm_compass(position = c(.58,.47))
#tmap_save(hotspots_max_mcf_peak, here("climate_change_data", "output_maps","hotspots_max_mcf_peak.png"))
subcatch_max_peak <- subcatch_max %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
ggplot() +
geom_sf(aes(fill = peak_runoff_cfs)) +
theme_bw() +
scale_fill_gradient(low = "#D1F2EB", high = "#148F77", name = "Peak Discharge (cfs)",
breaks = c(0, 25, 50, 75, 100 )) +
labs(title = "Max MCF") +
geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
# geom_sf_label(aes(label = subcatchment), size = 1.5) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90))+
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank(),
panel.border=element_blank()) +
theme_void() # for faculty presentation, remove lines and axis elements
subcatch_max_peak
# ggsave('subcatch_max_peak_removedsubs_themevoid.png', subcatch_max_peak, width = 5, height = 7, units = "in")
subcatch_max_total <- subcatch_max %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
ggplot() +
geom_sf(aes(fill = runoff_coeff)) +
theme_bw() +
scale_fill_gradient(low = "#D6EAF8", high = "#2874A6", name = "Runoff Coefficient",
breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1.0),
limits = c(0,1.0)) +
labs(title = "Max MCF") +
geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
# geom_sf_label(aes(label = subcatchment), size = 1.5) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank(),
panel.border=element_blank()) +
theme_void()
subcatch_max_total
# ggsave('subcatch_max_total.png', subcatch_max_total, width = 5, height = 7, units = "in")
Maps for MCF min storm hotspots
results_min_mcf<- read_csv("results_min_mcf.csv") ##read in file from code above
#Combine subcatchments outline with wet storm results
subcatch_min <- read_sf(dsn = here("climate_change_data", "SWMM_results","shapefiles"), layer = "subcatch_outline") %>%
st_transform(st_crs(4326)) %>%
clean_names() %>%
select(subcatchment = objectid_1) %>%
merge(results_min_mcf) %>%
filter(subcatchment != "5")
# Total volume hotspots
hotspots_min_mcf_total <- tm_basemap("OpenStreetMap.Mapnik") +
tm_shape(subcatch_min, unit = "Miles") +
tm_polygons("runoff_coeff", alpha = 0.8, palette = "Blues", style = "cont", n=8,
legend.hist = TRUE, title = "Runoff Coefficient") +
tm_layout(title = "March 2009 Min MCF storm", inner.margins=c(.05, .05, 0.1, .53),
legend.position = c(.6,.32), legend.title.size = 1.4, legend.text.size = 1) +
tm_text("subcatchment", size = 0.3) +
tm_scale_bar(position = c(.6,.59), breaks = c(0, 0.2, 0.4, 0.6, 0.8,1)) +
tm_compass(position = c(.58,.52))
##tmap_save(hotspots_wet_total, here("5.Results_Maps", "output_maps","hotspots_wet_total.png"))
# Peak flow hotspots
hotspots_min_mcf_peak <- tm_basemap("OpenStreetMap.Mapnik") +
tm_shape(subcatch_min, unit = "Miles") +
tm_polygons("peak_runoff_cfs", alpha = 0.75, palette = "Greens", style = "cont", n=8,
legend.hist = TRUE, title = "Peak Discharge (cfs)") +
tm_layout(title = "March 2009 Min MCF storm", inner.margins=c(.05, .05, 0.1, .53),
legend.position = c(.6,.27), legend.title.size = 1.4, legend.text.size = 1) +
tm_text("subcatchment", size = 0.3) +
tm_scale_bar(position = c(.6,.54), breaks = c(0, 0.2, 0.4, 0.6, 0.8,1)) +
tm_compass(position = c(.58,.47))
##tmap_save(hotspots_wet_peak, here("5.Results_Maps", "output_maps","hotspots_wet_peak.png"))
subcatch_min_peak <- subcatch_min %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
ggplot() +
geom_sf(aes(fill = peak_runoff_cfs)) +
theme_classic()+
scale_fill_gradient(low = "#D1F2EB", high = "#148F77", name = "Peak Discharge (cfs)",
breaks = c(0, 3, 6, 9, 12, 15)) +
labs(title = "Min MCF") +
geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
# geom_sf_label(aes(label = subcatchment), size = 1.5) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank(),
panel.border=element_blank()) +
theme_void() # for facily presentation, remove lines and axis elements
subcatch_min_peak
# ggsave('subcatch_min_peak_removedsub_themevoid.png', subcatch_min_peak, width = 5, height = 7, units = "in")
subcatch_min_total <- subcatch_min %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
ggplot() +
geom_sf(aes(fill = runoff_coeff)) +
theme_bw() +
scale_fill_gradient(low = "#D6EAF8", high = "#2874A6", name = "Runoff Coefficient",
breaks = c(0, 0.2, 0.4, 0.6, 0.8, 1.0),
limits = c(0,1.0)) +
labs(title = "Min MCF") +
geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
# geom_sf_label(aes(label = subcatchment), size = 1.5) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank(),
panel.border=element_blank()) +
theme_void()
subcatch_min_total
# ggsave('subcatch_min_total.png', subcatch_min_total, width = 5, height = 7, units = "in")
Normalized Peak Discharge Maps
subcatch_max_normalized <- read_sf(dsn = here("climate_change_data", "SWMM_results","shapefiles"), layer = "subcatch_outline") %>%
st_transform(st_crs(4326)) %>%
clean_names() %>%
select(subcatchment = objectid_1) %>%
merge(results_max_mcf_normalized) %>%
filter(subcatchment != "5")
## normalized Max MCF Peak Discharge
subcatch_max_peak_normalized <- subcatch_max_normalized %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
ggplot() +
geom_sf(aes(fill = peak_runoff_cfs_norm)) +
theme_bw() +
scale_fill_gradient(low = "#D1F2EB", high = "#148F77", name = "Normalized Peak Discharge (cfs/sqft)") +
labs(title = "Max MCF") +
geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
# geom_sf_label(aes(label = subcatchment), size = 1.5) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90))+
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank(),
panel.border=element_blank()) +
theme_void() # for faculty presentation, remove lines and axis elements
subcatch_max_peak_normalized
# ggsave('subcatch_max_peak_removedsubs_themevoid_norm.png', subcatch_max_peak_normalized, width = 5, height = 7, units = "in")
subcatch_min_normalized <- read_sf(dsn = here("climate_change_data", "SWMM_results","shapefiles"), layer = "subcatch_outline") %>%
st_transform(st_crs(4326)) %>%
clean_names() %>%
select(subcatchment = objectid_1) %>%
merge(results_min_mcf_normalized) %>%
filter(subcatchment != "5")
## normalized Min MCF Peak Discharge
subcatch_min_peak_normalized <- subcatch_min_normalized %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
ggplot() +
geom_sf(aes(fill = peak_runoff_cfs_norm)) +
theme_classic()+
scale_fill_gradient(low = "#D1F2EB", high = "#148F77", name = "Normalized Peak Discharge (cfs/sqft)") +
labs(title = "Min MCF") +
geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 1.8) +
# geom_sf_label(aes(label = subcatchment), size = 1.5) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank(),
panel.border=element_blank()) +
theme_void() # for facily presentation, remove lines and axis elements
subcatch_min_peak_normalized
# ggsave('subcatch_min_peak_removedsub_themevoid_norm.png', subcatch_min_peak_normalized, width = 5, height = 7, units = "in")
binned runoff coefficient maps for min and max MCF
library(RColorBrewer)
library(metR)
# View a single RColorBrewer palette by specifying its name
display.brewer.pal(n = 7, name = 'Blues')
# Hexadecimal color specification
brewer.pal(n = 7, name = "Blues")
## [1] "#EFF3FF" "#C6DBEF" "#9ECAE1" "#6BAED6" "#4292C6" "#2171B5" "#084594"
mycols <- brewer.pal(n = 11, name = "Blues")
mybreaks_max<- c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
subcatch_min_total <- subcatch_min %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
mutate(runoff_coeff_cut = cut_number(runoff_coeff, n = 10)) %>%
ggplot() +
geom_sf(aes(fill = runoff_coeff_cut)) +
theme_bw() +
scale_fill_gradientn(colours = mycols,
breaks= mybreaks_max,
super = metR::ScaleDiscretised,
name = "Runoff Coefficient") +
theme(legend.position = "bottom") +
labs(title = "Min MCF") +
geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 3.1) +
# geom_sf_label(aes(label = subcatchment), size = 1.5) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank(),
panel.border=element_blank()) +
theme_void()
subcatch_min_total
# ggsave('subcatch_min_total_binned.png', subcatch_min_total, width = 5, height = 7, units = "in")
mycols <- brewer.pal(n = 11, name = "Blues")
mybreaks_max<- c(0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9)
subcatch_max_total <- subcatch_max %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
mutate(runoff_coeff_cut = cut_number(runoff_coeff, n = 10)) %>%
ggplot() +
geom_sf(aes(fill = runoff_coeff_cut)) +
theme_bw() +
scale_fill_gradientn(colours = mycols,
breaks= mybreaks_max,
super = metR::ScaleDiscretised,
name = "Runoff Coefficient") +
labs(title = "Max MCF") +
geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 3.1) +
# geom_sf_label(aes(label = subcatchment), size = 1.5) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank(),
panel.border=element_blank()) +
theme_void()
subcatch_max_total
# ggsave('subcatch_max_total_binned.png', subcatch_max_total, width = 5, height = 7, units = "in")
binned peak flow maps for min and max MCF
mycols2 <- brewer.pal(n = 11, name = "Greens")
mybreaks2_min<- c(0.0, 2.0, 4.0, 6.0, 8.0, 10.0, 12.0, 14.0, 16.0)
subcatch_min_peak <- subcatch_min %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
mutate(peak_runoff_cfs_cut = cut_number(peak_runoff_cfs, n = 11)) %>%
ggplot() +
geom_sf(aes(fill = peak_runoff_cfs_cut)) +
theme_classic()+
scale_fill_gradientn(colours = mycols2,
breaks= mybreaks2_min,
super = metR::ScaleDiscretised,
name = "Peak Discharge (cfs)") +
labs(title = "Min MCF") +
geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 3.1) +
# geom_sf_label(aes(label = subcatchment), size = 1.5) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90)) +
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank(),
panel.border=element_blank()) +
theme_void() # for facily presentation, remove lines and axis elements
subcatch_min_peak
# ggsave('subcatch_min_peak_removedsub_themevoid_binned.png', subcatch_min_peak, width = 5, height = 7, units = "in")
mycols3 <- brewer.pal(n = 11, name = "Greens")
mybreaks3_max<- c(0.0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)
subcatch_max_peak <- subcatch_max %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
mutate(peak_runoff_cfs_cut = cut_number(peak_runoff_cfs, n = 11)) %>%
ggplot() +
geom_sf(aes(fill = peak_runoff_cfs_cut)) +
theme_classic()+
scale_fill_gradientn(colours = mycols3,
breaks= mybreaks3_max,
super = metR::ScaleDiscretised,
name = "Peak Discharge (cfs)") +
labs(title = "Max MCF") +
geom_sf_text(aes(label = subcatchment), colour = "black", fontface = "bold", size = 3.1) +
# geom_sf_label(aes(label = subcatchment), size = 1.5) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
theme(axis.text.x = element_text(angle = 90))+
theme(axis.text.x = element_blank(),
axis.text.y = element_blank(),
axis.ticks = element_blank(),
rect = element_blank(),
panel.border=element_blank()) +
theme_void() # for faculty presentation, remove lines and axis elements
subcatch_max_peak
# ggsave('subcatch_max_peak_removedsubs_themevoid_binned.png', subcatch_max_peak, width = 5, height = 7, units = "in")
combine max and min runoff coefficient results into single table showing the catchment #, min, and max runoff coefficient, slope, percent impervious, and area (sqft)
results_min_mcf_rc <- results_min_mcf %>%
select(subcatchment, Slope, percent_imp, Area_sqft, runoff_coeff) %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
arrange(desc(runoff_coeff)) %>%
top_n(20,runoff_coeff)
results_max_mcf_rc <- results_max_mcf %>%
select(subcatchment, Slope, percent_imp, Area_sqft, runoff_coeff) %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
arrange(desc(runoff_coeff)) %>%
top_n(20,runoff_coeff)
results_runoffcoef_table <- results_min_mcf_rc %>%
full_join(results_max_mcf_rc, by=c("subcatchment", "Slope", "percent_imp", "Area_sqft")) %>%
kable(col.names=c("Subcatchment","Slope", "Percent Impervious", "Area (sqft)", "Runoff Coefficient, Min MCF",
"Runoff Coefficient, Max MCF")) %>%
kable_styling(bootstrap_options = "striped")
results_runoffcoef_table
|
Subcatchment
|
Slope
|
Percent Impervious
|
Area (sqft)
|
Runoff Coefficient, Min MCF
|
Runoff Coefficient, Max MCF
|
|
45
|
10.837372
|
75.81188
|
205273.48
|
0.594
|
0.815
|
|
47
|
9.017865
|
73.67874
|
163180.85
|
0.586
|
0.799
|
|
60
|
10.200056
|
66.27512
|
130703.86
|
0.531
|
0.821
|
|
40
|
7.774446
|
65.26444
|
78272.63
|
0.528
|
0.762
|
|
51
|
3.205696
|
65.70618
|
194382.01
|
0.516
|
0.746
|
|
54
|
11.985041
|
65.90850
|
301206.44
|
0.514
|
0.748
|
|
21
|
2.133592
|
64.38962
|
78832.74
|
0.511
|
NA
|
|
29
|
4.303893
|
63.75090
|
62214.52
|
0.510
|
0.747
|
|
67
|
13.094514
|
65.31544
|
376961.02
|
0.509
|
0.809
|
|
71
|
1.229756
|
63.33945
|
161251.17
|
0.493
|
0.793
|
|
5
|
4.973908
|
59.22529
|
42008.07
|
0.489
|
0.809
|
|
49
|
11.626618
|
61.33458
|
317060.05
|
0.484
|
NA
|
|
46
|
9.180995
|
59.97949
|
209977.29
|
0.481
|
0.789
|
|
13
|
2.801558
|
60.72206
|
210147.80
|
0.473
|
NA
|
|
59
|
9.355130
|
58.61615
|
241399.55
|
0.465
|
0.751
|
|
22
|
5.224538
|
57.36459
|
106889.43
|
0.463
|
0.745
|
|
35
|
2.502320
|
58.10591
|
163977.29
|
0.461
|
NA
|
|
36
|
3.413503
|
57.37835
|
103198.24
|
0.458
|
NA
|
|
63
|
11.350848
|
59.59508
|
1260905.98
|
0.454
|
0.746
|
|
23
|
2.091764
|
60.15849
|
712111.05
|
0.452
|
0.759
|
|
68
|
16.778001
|
56.40240
|
306366.86
|
NA
|
0.813
|
|
65
|
12.285943
|
54.98557
|
309782.87
|
NA
|
0.770
|
|
14
|
10.062202
|
51.97367
|
244581.39
|
NA
|
0.764
|
|
91
|
7.225773
|
49.19777
|
543493.09
|
NA
|
0.755
|
|
89
|
12.267217
|
54.37988
|
929708.78
|
NA
|
0.748
|
combine max and min peak runoff (cfs) results into single table showing the catchment #, min, and max peak runoff (cfs)
results_min_mcf_peak <- results_min_mcf %>%
select(subcatchment, Slope, percent_imp, Area_sqft, peak_runoff_cfs) %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
arrange(desc(peak_runoff_cfs)) %>%
top_n(20,peak_runoff_cfs)
results_max_mcf_peak <- results_max_mcf %>%
select(subcatchment, Slope, percent_imp, Area_sqft, peak_runoff_cfs) %>%
filter(!subcatchment %in% c(1, 2,94,3,38, 24)) %>%
arrange(desc(peak_runoff_cfs)) %>%
top_n(20,peak_runoff_cfs)
results_peakflow_table <- results_min_mcf_peak %>%
full_join(results_max_mcf_peak, by=c("subcatchment", "Slope", "percent_imp", "Area_sqft")) %>%
kable(col.names=c("Subcatchment","Slope", "Percent Impervious", "Area (sqft)","Peak Runoff (cfs), Min MCF",
"Peak Runoff (cfs), Max MCF")) %>%
kable_styling(bootstrap_options = "striped")
results_peakflow_table
|
Subcatchment
|
Slope
|
Percent Impervious
|
Area (sqft)
|
Peak Runoff (cfs), Min MCF
|
Peak Runoff (cfs), Max MCF
|
|
63
|
11.3508475
|
59.595082
|
1260906.0
|
15.29
|
106.39
|
|
28
|
9.6951530
|
45.648739
|
1400355.0
|
13.56
|
99.65
|
|
42
|
12.7741916
|
52.075244
|
1097149.0
|
12.15
|
80.56
|
|
11
|
0.9475857
|
60.176841
|
1461121.3
|
11.91
|
105.54
|
|
89
|
12.2672167
|
54.379880
|
929708.8
|
10.85
|
79.13
|
|
74
|
0.5130342
|
49.951879
|
1593473.1
|
10.07
|
90.09
|
|
62
|
4.9693069
|
57.621143
|
816823.7
|
9.66
|
62.88
|
|
58
|
13.4737125
|
50.489651
|
794612.3
|
8.76
|
57.25
|
|
78
|
0.6625196
|
45.869531
|
1272102.3
|
8.38
|
69.55
|
|
23
|
2.0917640
|
60.158489
|
712111.1
|
8.09
|
59.82
|
|
76
|
12.6689327
|
21.946154
|
1533842.8
|
7.40
|
52.96
|
|
7
|
24.3708501
|
14.101788
|
2174275.2
|
6.79
|
56.98
|
|
75
|
3.1037990
|
49.969260
|
646779.6
|
6.65
|
47.57
|
|
79
|
0.7809627
|
55.674414
|
703208.6
|
6.63
|
48.58
|
|
91
|
7.2257735
|
49.197766
|
543493.1
|
5.86
|
47.33
|
|
52
|
14.2092259
|
49.068201
|
521285.7
|
5.66
|
38.27
|
|
10
|
13.0349991
|
46.827316
|
527259.5
|
5.46
|
38.42
|
|
67
|
13.0945141
|
65.315445
|
376961.0
|
5.35
|
37.74
|
|
34
|
6.5716791
|
55.845327
|
439337.8
|
5.31
|
35.34
|
|
55
|
14.0124563
|
43.238717
|
456026.8
|
4.37
|
NA
|
|
20
|
28.7382083
|
5.926749
|
3260351.4
|
NA
|
64.99
|
comparing historical vs Max MCF runoff coeffiecients from Dornan et al., 2020
## runoff coefficients differences
dornan_results_top20_rc <- read_csv("~/Desktop/Aloha-Aina-Master/climate_change_data/SWMM_results/dornan_results_top20_rc.csv") %>%
select(subcatchment, runoff_coeff) # read in dorana et al RC results
results_max_mcf_rc_summary <- results_max_mcf_rc %>% select(subcatchment, runoff_coeff) # select just our subcatchment and runoff coefficients to be able to subtract
combined_rc_results <- dornan_results_top20_rc %>% right_join(results_max_mcf_rc_summary, by = c("subcatchment")) # combine dornan results with our max MCF reslts
combined_rc_results = combined_rc_results %>% mutate(subtract = runoff_coeff.y-runoff_coeff.x) # add column with difference between max MCF rc results - dornan historical rc results with March 14, 09 storm
## peak flow differences
dornan_results_top20_peakflow <- read_csv("~/Desktop/Aloha-Aina-Master/climate_change_data/SWMM_results/dornan_results_top20_peakflow.csv") %>%
select(subcatchment, peak_runoff_cfs) # read in Team Kahuawais peak flow results
results_max_mcf_peak_summary <- results_max_mcf_peak %>% select(subcatchment, peak_runoff_cfs) # select just our subcatchment and peak runoff cfs to be able to subtract
combined_peak_results <- dornan_results_top20_peakflow %>% right_join(results_max_mcf_peak_summary, by = c("subcatchment"))
combined_peak_results = combined_peak_results %>% mutate(subtract = peak_runoff_cfs.y-peak_runoff_cfs.x)